home *** CD-ROM | disk | FTP | other *** search
/ The X-Philes (2nd Revision) / The X-Philes Number 1 (1995).iso / xphiles / hp48_1 / fft_jm < prev    next >
Internet Message Format  |  1995-03-31  |  12KB
  1. From: John Paul Morrison <jmorriso@ee.ubc.ca>
  2. Subject:  v03i043:  fft_jm - New FFT v1.0, Part01/01
  3. Newsgroups: comp.sources.hp48
  4. Followup-To: comp.sys.hp48
  5. Approved: spell@seq.uncwil.edu
  6.  
  7. Checksum: 3095338785 (verify with brik -cv)
  8. Submitted-by: John Paul Morrison <jmorriso@ee.ubc.ca>
  9. Posting-number: Volume 3, Issue 43
  10. Archive-name: fft_jm/part01
  11.  
  12.  
  13. BEGIN_DOC fft.doc
  14. (main stuff: FFT, IFFT, the rest are just token utilities)
  15. RPL FFT program
  16.  
  17. This is an FFT program written (mostly) in RPL. I've included
  18. some support programs that will help analysis of signals etc.
  19.  
  20. The original algorithm and code (VFFT) was written by Nicola Catacchio
  21. (e-mail: ares@alessia.unipd.it). Her code was the fastest user code
  22. I have seen, so that's why I recoded it in RPL. This new FFT
  23. code is 1.7 time faster for complex vectors, and 1.5 times faster
  24. for real vectors, than the original VFFT program, and at least twice
  25. as fast as others.
  26.  
  27. some sample times:
  28.  
  29. 128 element (complex)
  30. RPL FFT: 49 s
  31. VFFT:    84 s
  32. others: 100 s
  33.  
  34. 64 element (complex)
  35. RPL FFT: 20 s
  36. VFFT:    35 s
  37. others:  41 s
  38.  
  39. RPL code is often around 10 times faster than user code, but
  40. unfortunately, the bottle-neck here is with the RPL complex multiply,
  41. and less significantly with other complex operations. The HP routine
  42. for a complex multiply is done in RPL, not machine code (ie they call
  43. the rpl multiply, which is in machine code).
  44.  
  45. To get some more performance, I recoded the complex mulitplies, since
  46. the built in complex multiply uses extended precision arithmetic.
  47. Standard precision is faster, and the error is not significant (since
  48. the FFT algorithm doesn't compund errors the way other algorithms do).
  49.  
  50. If anyone has looked at the FFT problem, I'd like to talk to them,
  51. and try and improve performance. Send e-mail to John Paul Morrison,
  52. jmorriso@ee.ubc.ca. I have not posted the RPL source-code, but I
  53. can email it to those interested.
  54.  
  55. Rather than waiting any longer (until the routines are elegantly coded)
  56. here are the programs, with brief descriptions:
  57.  
  58. FFT
  59.     This takes a real or complex vector. If the dimensions are not
  60.     2^N it will give an error message.
  61.  
  62.     (run ASC\-> to decode after downloading)
  63.  
  64. IFFT
  65.     Inverts the FFT.
  66.  
  67. SAMPLE
  68.     'algebraic'
  69.     'name'
  70.     start
  71.     finish
  72.     n
  73.  
  74.     This returns n samples in a vector. It samples the 'algebraic' with
  75.     independent variable 'name' from times start to finish.
  76.  
  77. PLOT
  78.     Useful for time domain plots.
  79.     Does a bar-plot of a vector V(1...N). V(1) is plotted on the left,
  80.     and V(N) is on the right.
  81.  
  82. PLOTF
  83.     Useful for frequency domain plots.
  84.     Does a bar-plot of a vector V(1...N). V(1) is in the center, plotted
  85.     up to V(N/2) on the right. V(N/2+1)..V(N). are on the left.
  86.  
  87. CLR
  88.     Clears the pict, and statistics information etc.
  89.  
  90. SINC
  91. RECT
  92. TRI
  93.     Defines some familiar communications signals.
  94.  
  95. LIN
  96.     LIN(t,t1,f1,t2,f2)
  97.     This defines a straight line from (t1,f1) to (t2,f2) within the range
  98.     t1...t2; with zero otherwise. This is useful for defining discontinuous
  99.     functions with lots of straight lines in it.
  100.  
  101.  
  102. John Paul Morrison.
  103. END_DOC
  104.  
  105. BEGIN_RDME fft.rdm
  106.     [ The rpl version of this program has the program fft asc'ed.
  107.     Fft has to be converted before it can be used.  However, the
  108.     asc & uuencoded parts of this post already have fft
  109.     converted so you do not have to convert it yourself.
  110.  
  111.     Chris ]
  112. END_RDME
  113.  
  114. BYTES: #FF8Ch 2003
  115.  
  116. BEGIN_RPL fft.rpl
  117. %%HP: T(3)A(R)F(.);
  118. DIR
  119.   PLOT
  120.     \<< DUP OBJ\-> EVAL
  121. \-> N
  122.       \<< { N 1 }
  123. \->ARRY STO\GS BARPLOT
  124. GRAPH
  125.       \>>
  126.     \>>
  127.   PLOTF
  128.     \<< DUP OBJ\-> EVAL
  129. \-> N
  130.       \<< 1 N 2 /
  131.         START N
  132. ROLL
  133.         NEXT { N 1
  134. } \->ARRY STO\GS
  135. BARPLOT GRAPH
  136.       \>>
  137.     \>>
  138.   ABSV
  139.     \<< OBJ\-> EVAL \-> N
  140.       \<< 1 N
  141.         START N
  142. ROLL ABS
  143.         NEXT N
  144. \->ARRY
  145.       \>>
  146.     \>>
  147.   SAMPLE
  148.     \<< 0 \-> f t t0 t1
  149. n s
  150.       \<< { } t + 'f'
  151. APPLY f = DEFINE t1
  152. t0 - n 1 - / 's'
  153. STO 1 n
  154.         FOR i t0 f
  155. EVAL s 't0' STO+
  156.         NEXT n
  157. \->ARRY
  158.       \>>
  159.     \>>
  160.   SINC
  161.     \<< \-> t
  162.       \<<
  163.         IF t 0 \=/
  164.         THEN t \pi *
  165. DUP SIN SWAP /
  166.         ELSE 1
  167.         END
  168.       \>>
  169.     \>>
  170.   TRI
  171.     \<< \-> t
  172.       \<< t ABS
  173.         IF DUP 1 \<=
  174.         THEN 1 -
  175. NEG
  176.         ELSE DROP 0
  177.         END
  178.       \>>
  179.     \>>
  180.   RECT
  181.     \<< \-> t
  182.       \<< t ABS .5 \<=
  183.       \>>
  184.     \>>
  185.   LIN
  186.     \<< \-> t t0 f0 t1
  187. f1
  188.       \<<
  189.         IF t t0 < t
  190. t1 > OR
  191.         THEN 0
  192.         ELSE f0 t
  193. t1 - * f1 t t0 - *
  194. - t0 t1 - /
  195.         END
  196.       \>>
  197.     \>>
  198.   FFT
  199. "D9D20E16327792000000000000000100000000000000000EEDA1290D19C2A26C
  200. 7D178BF1F49B1ED2A2F49B150FA13CE223ABB14B2A2D9AE1AFE22D9D20900D1E
  201. 4A20510001050000000000000933A1B21305BF22D9D20AEC81881303004038D3
  202. 0CB916D9D2088130E8E3088130FC436FED30D00405923062E267F37088130122
  203. 70E0E3038D30CB916D9D2012270523302C230353161F03612270E9330B21306B
  204. 316322302A170B9826EE170D9D202C23021E26E8E3032230B21305E170CBD304
  205. 337079470B21309FF302C2309034647A2003D4303D4303D4303D4303D43B2130
  206. 0D470A69159FF3059230FBD81E6BA2ED2A2E6BA2EF9A2AEC8162E267F3707F42
  207. 583416D2E3088130C54160ED30FED300E5160F5162A1707E316006162A170E04
  208. 16C5416CC2162C230C2D5059230C2D50D7C26CB9A24C016EF116CB9A2E9016EF
  209. 116CB9A2CAF06CB9A2BBF06189A2CAF06479A272C50E041654D26CA130CFC15C
  210. AF0661C15E0416E9330E0416C5416C2316E0416C5416F6E308F7268813000616
  211. D00404EC3032230119200000838D3057B308C1702C230C2D5059230C2D50D7C2
  212. 6CB9A24C016EF116CB9A2E9016EF116CB9A2CAF06CB9A2BBF06189A2CAF06479
  213. A272C507E316A18260F5166B3164EC308C17044230C5416F6E30526162BB1570
  214. 1252BB1543370442308341679470FBD81900D1B21305DF2293632B21304AA0"
  215.   IFFT
  216.     \<< CONJ FFT CONJ
  217. DUP SIZE EVAL /
  218.     \>>
  219.   CLR
  220.     \<< { \GSPAR PPAR
  221. \GSDAT PICT } PURGE
  222.     \>>
  223.   CST { PLOT PLOTF
  224. ABSV SAMPLE SINC
  225. TRI RECT LIN FFT
  226. IFFT C\->R CLR }
  227. END
  228. END_RPL
  229.  
  230. BEGIN_ASC fft.asc
  231. %%HP: T(3)A(D)F(.);
  232. "69A20FF72FA0000000303435453047A2084E204005C4F44584E205005C4F4456
  233. 484E20401424356584E20603514D405C45484E20403594E43484E20304525948
  234. 4E20402554344584E2030C494E484E203064644584E204094646445E89C184E2
  235. 03034C425B2130CB0003034C42530D9D20E163247A2084E20405805142584E20
  236. 400505142584E204058441445634E1B2130EFE0293632B2130F5000409464644
  237. 540D9D20E1632E6AA184E2030646445E6AA178BF18B9C1EB3A150FA193632B21
  238. 30B40003064644530D9D20E16327792000000000000000100000000000000000
  239. EEDA1290D19C2A26C7D178BF1F49B1ED2A2F49B150FA13CE223ABB14B2A2D9AE
  240. 1AFE22D9D20900D1E4A20510001050000000000000933A1B21305BF22D9D20AE
  241. C81881303004038D30CB916D9D2088130E8E3088130FC436FED30D0040592306
  242. 2E267F3708813012270E0E3038D30CB916D9D2012270523302C230353161F036
  243. 12270E9330B21306B316322302A170B9826EE170D9D202C23021E26E8E303223
  244. 0B21305E170CBD304337079470B21309FF302C2309034647A2003D4303D4303D
  245. 4303D4303D43B21300D470A69159FF3059230FBD81E6BA2ED2A2E6BA2EF9A2AE
  246. C8162E267F3707F42583416D2E3088130C54160ED30FED300E5160F5162A1707
  247. E316006162A170E0416C5416CC2162C230C2D5059230C2D50D7C26CB9A24C016
  248. EF116CB9A2E9016EF116CB9A2CAF06CB9A2BBF06189A2CAF06479A272C50E041
  249. 654D26CA130CFC15CAF0661C15E0416E9330E0416C5416C2316E0416C5416F6E
  250. 308F7268813000616D00404EC3032230119200000838D3057B308C1702C230C2
  251. D5059230C2D50D7C26CB9A24C016EF116CB9A2E9016EF116CB9A2CAF06CB9A2B
  252. BF06189A2CAF06479A272C507E316A18260F5166B3164EC308C17044230C5416
  253. F6E30526162BB15701252BB1543370442308341679470FBD81900D1B21305DF2
  254. 293632B21304040030C494E430D9D20E16321C432D6E201047D6E20204703D6E
  255. 20206603D6E20204713D6E20206613E16323CE22D6E201047D6E20204703EBBE
  256. 1D6E201047D6E20204713D5CE1908E1AFE224B2A25BF22D9D20D6E20206603D6
  257. E201047D6E2020471390DA1EEDA1D6E20206613D6E201047D6E2020470390DA1
  258. EEDA190DA1D6E20204703D6E2020471390DA150FA1B21305DF22EF53293632B2
  259. 13033100402554344540D9D20E16321C432D6E201047E1632D6E201047F1AA13
  260. 39209990000000000050CFCE1EF53293632B2130060003045259430D9D20E163
  261. 21C432D6E201047E1632D6E201047F1AA13CE2278BF19C2A2CFCE1AFE22D9D20
  262. 9C2A290DA1599A1B21305BF22D9D208DBF14B2A2B21305DF22EF53293632B213
  263. 049000403594E43440D9D20E16321C432D6E201047E16323CE22D6E2010474B2
  264. A2D9AE1AFE22D9D20D6E201047DBAA1EEDA178BF1CA4B1DBBF150FA1B21305BF
  265. 229C2A25DF22EF53293632B213059000603514D405C45460D9D20E16324B2A21
  266. C432D6E201066D6E201047D6E20204703D6E20204713D6E2010E6D6E201037E1
  267. 63247A20B2130D6E20104776BA145632D6E20106697632D55F1D6E2010668D8A
  268. 156D02D6E20204713D6E2020470390DA1D6E2010E69C2A290DA150FA145632D6
  269. E20103797632DCC029C2A2D6E2010E60A132D6E201096D6E20204703D6E20106
  270. 6EB3A1D6E20103745632D6E2020470397632B4402C4232D6E2010E6900D1EF53
  271. 293632B2130B6100401424356540D9D20E1632B7FC1EB3A11C432D6E2010E4E1
  272. 6329C2A2D6E2010E430132D6E2010E45BCF1F1AA1C4232D6E2010E4900D1EF53
  273. 293632B2130B70005005C4F4456450D9D20E163278BF1B7FC1EB3A11C432D6E2
  274. 010E4E16329C2A2D6E2010E4ED2A250FA130132D6E2010E45BCF1C423247A20D
  275. 6E2010E49C2A2B2130900D1B0DF133102AB2E1EF53293632B21305A0004005C4
  276. F44540D9D20E163278BF1B7FC1EB3A11C432D6E2010E4E163247A20D6E2010E4
  277. 9C2A2B2130900D1B0DF133102AB2E1EF53293632B2130310F"
  278. END_ASC
  279.  
  280. BYTES: #F013h 1470
  281.  
  282. BEGIN_UU fft.uue
  283. begin 644 fft
  284. M2%!(4#0X+466*O!_\@H````#0U-4`W0J@.0"!%!,3U1(+E``Q?1$983D`@1!$
  285. M0E-62"Y@,!74!,54A.0"!%-)3D-(+C!`)96$Y`($4D5#5$@N,,"4Y(3D`@-&<
  286. M1E1(+D"09&1$Y9@<2"XP,,0DM1(#O``P,,0D-=#9`AXV0J<"2"Y`4`@5)(7D=
  287. M`@104$%22"Y`4$@41&5#'BLQX.\@.3:R$@-?`$"09&1$1=#9`AXVXJ8:2"XP0
  288. M8&1$Y:8:A_N!FQR^HU'P&CDVLA(#2P`P8&1$-=#9`AXV<I<"``````````$`>
  289. M`````````.ZM(0D=R:)B?!V'^_&4&]ZB\I0;!:\Q["*CNT$K*IWJH>\BG2V0B
  290. M`!U.*E`!``$%````````.:.Q$@.U+]+9`NJ,@1@#`T`PV`.\&=;9`H@QX.@#,
  291. MB#'P3&/O/=``!)4R8.)B]W.`&`,A<N#@`X,]P)MAG2T0(@<E,R`L`U,3%@]C3
  292. M(7+@.0,K,6`[82,R(!H'FRCF'@>=+2`L`Q(NYN@#(S*P$@/E<<#;`S1S<$D'^
  293. M*S&0_P/",I`P9'0J`-,T,$T#TS0P30/3-"LQ`$T':AF5_P.5,O#;&&ZKXBTJ'
  294. M;JOBGRKJC&'B8O=S<$]2.!36X@.(,<!%8>`]\-X#X!4&7V&B<7`^80`6)AH'W
  295. M#A3&16',$B8L`RQ=4"D#+%W0QV*\J4(,8?X1QILJGA#F'V&\J<+Z8+RILOM@%
  296. M@:G"^F!TJ7+"!0X45M1BK#'`SU&L#V;!40X4YCD##A3&16$L$^9`85P4]N8#9
  297. M^">&&`,`%M8`!.0\,"(#$2D``("#/5"W`\AQ("P#+%U0*0,L7=#'8KRI0@QAH
  298. M_A'&FRJ>$.8?8;RIPOI@O*FR^V"!J<+Z8'2I<L(%YQ.F@6+P%68[8>0\@!P'E
  299. M1#+`16%O/E!B8;(;=1!2LAM%,P=$,H!#89=T\-L8"="Q$@/5+Y)C(RLQ0$```
  300. M`TQ)3@.=+>!A(\$TTN8"`71M+B!`!]/F`@)F,&TN($`7T^8"`F8Q'C8R[")M-
  301. M+A!`U^8"`G0PONO1Y@(!=&TN($`7T\4>">BA[R*THE+[(ITMT.8"`F8P;2X0<
  302. M0-?F`@)T,0FMX=X:;2X@8!;3Y@(!=&TN($`'D]`:[JV1T!IM+B!`!]/F`@)T+
  303. M,0FM4?`:*S%0_2+^-9)C(RLQ,!,`!%)%0U0$G2W@82/!--+F`@%T'C;2Y@(!^
  304. M=!^J,9,"F0D```````7\[.%?(SDVLA(#8``P0"65--#9`AXV$DPC;2X00.=A6
  305. M(VTN$$#WH1K#+G*X'\FBPL\>^B[2V0+)HI+0&I6IL1(#M2_2V0+8^T$K*BLQ*
  306. M4/TB_C628R,K,4`)``1324Y#!)TMX&$CP332Y@(!=!XV,NPB;2X00$<K*IWJ)
  307. MH>\BG2W0Y@(!=+VJX=X:A_O!2AN]^U'P&BLQ4/LBR:)2_2+^-9)C(RLQ4`D`E
  308. M!E-!35!,10:=+>!A([2B$DPC;2X08-;F`@%T;2X@0`?3Y@("=#%M+A#@UN8"Z
  309. M`7,>-D*G`BLQT.8"`71GJT%E(VTN$&"69R-=]='F`@%FV*A1UB!M+B!`%]/FZ
  310. M`@)T,`FMT>8"`6[)HI+0&@6O064C;2X0,)=G(\T,DBPJ;2X0X`8:(VTN$)#6-
  311. MY@("=#!M+A!@YCL:;2X0,$=E(VTN($`'DV<C2P3")"-M+A#@E@`=_C628R,K+
  312. M,;`6``1!0E-6!)TMX&$C>\_A.QK!--+F`@%.'C:2+"IM+A#@-!`C;2X0X%3+C
  313. M'Q^JP20C;2X0X)0`'?XUDF,C*S&P!P`%4$Q/5$8%G2W@82.'^['W'+ZC$4PC^
  314. M;2X0X.1A(\FBTN8"`4[>HE+P&@,QTN8"`4ZU_,$D(W0JT.8"`4[)HK(2`PG0`
  315. ML=`?,P&B*Q[^-9)C(RLQ4`H`!%!,3U0$G2W@82.'^['W'+ZC$4PC;2X0X.1AA
  316. ?(W0JT.8"`4[)HK(2`PG0L=`?,P&B*Q[^-9)C(RLQ`"LQ)
  317. ``
  318. end
  319. END_UU
  320. -- 
  321. __________________________________________________________________________
  322.  John Paul Morrison                     |
  323.  University of British Columbia, Canada |
  324.  Electrical Engineering                 |   .sig closed for repairs 
  325.                                         |             
  326.  jmorriso@ee.ubc.ca                     |
  327. ________________________________________|_________________________________
  328.  
  329.